Nucleation of a sodium droplet on Ceo 
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We investigate theoretically the progressive coating of Ceo by several sodium atoms. Density 
functional calculations using a nonlocal functional are performed for NaCeo and Na2C6o in various 
configurations. These data are used to construct an empirical atomistic model in order to treat larger 
sizes in a statistical and dynamical context. Fluctuating charges are incorporated to account for 
charge transfer between sodium and carbon atoms. By performing systematic global optimization 
in the size range 1 < n < 30, we find that Na„C6o is homogeneously coated at small sizes, and that a 
growing droplet is formed above n > 8. The separate effects of single ionization and thermalization 
are also considered, as well as the changes due to a strong external electric field. The present results 
are discussed in the light of various experimental data. 

PACS numbers: 61.48.+c,34.70.+e,36.40.Qv 

I. INTRODUCTION 

Bulk compounds of Ceo with alkali or alkaline earth metals can develop particularly interesting properties such 
as superconductivity. 1-4 At the gas phase level, the interaction between a fullerene molecule and metal atoms has 
attracted a significant attention, from both theoreticians 5-16 and experimentalists. 17-29 Coverage of fullerenes with 
semiconductor atoms or clusters has also been investigated. 30,31 For alkaline earth materials it is now commonly 
accepted that many atoms homogeneously surround the Cgo molecule in spherically centred shells. For other metals, 
q . such as lithium 5,11,12,21 a strong charge transfer takes place with carbon and the resulting ionic interaction may be 
stronger than the metallic bonding. As a consequence, homogeneous coating of Cgo is also observed at low coverage. 
The same situation is seen for potassium 6,11,18,20 and rubidium 20,23 adsorbed on Ceo- The adsorption of gold leads to 
a somewhat different picture, where the metal cluster is formed next to the fullerene. 28 In this case metallic bonding 
dominates over the partially ionic-covalent Au-C interactions, 
■^j- ■ Here we will focus on sodium, which has been the subject of a debate. 18,26,28,29 The results of mass spectrometry 
| measurements, and the magic numbers inferred from these measurements by the group of Martin, 18 lead to the con- 
^^O ■ elusion of a regular coating of Cgo- Photoelectron spectroscopy measurements by Palpant et aL 26,28 were interpreted 
rather similarly, these authors also providing some evidence that coating proceeds by trimers rather than monomers. 
On the contrary, the measurement of electric dipoles and polarizabilities by Dugourd et al. 29 seems to indicate that a 
segregated metal droplet is formed on the surface of the fullerene. Since the experimental conditions for the production 
of these clusters are similar in the three apparatus, then mainly the interpretations differ. 

Up to now, there have been essentially three theoretical investigations of this problem, at least to our knowledge. 7,9,16 
The work by Rubio et al. 7 assumes a complete wetting of Cgo by sodium with a phenomenological two-shell jcllium 
I ■ description. This does not help much in elucidating the actual structure of Na„Ceo compounds. The ab initio 
calculations by Hira and Ray, 9 performed at the unrestricted Hartree-Fock (UHF) level, and the recent density 
functional theory (DFT) calculations by Hamamoto and coworkers, 16 using the local density approximation (LDA), 
predict quite different results for the physical properties and favored structures of Na2 on Cgo ■ 

In order to address this problem, we have constructed an empirical atomistic model, allowing extensive sampling of 
the configuration space in a wide range of sizes. In turn, this enhanced sampling enables one to achieve unconstrained 
global optimization as well as dynamical and finite temperature studies. 

For this model to be physically and chemically realistic, we first carefully performed first-principles calculations of 
the binding of Na and Na2 on Cgo- These calculations, described and discussed in the next section, are converted 
into a suitable set of parameters for our model, as fully described in Sec. III. The main results can be separated 
into five categories, which will be presented in the following order. After describing the most stable structures found 
with the present model, and some of their properties, we discuss the distinct effects of charging and thermalizing 
the clusters, and the possible influence of a strong external electric field. The diffusion dynamics of the metal atoms 
over the Cgo surface is also discussed. The main results are summarized and discussed in Sec. V, in the light of the 
various experimental and theoretical works available. A tentative rationalization for the apparently contradictory 
observations is proposed. 32 
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II. FIRST-PRINCIPLES CALCULATIONS 
A. Methods 

Density functional theory calculations have been used to provide reference data for the interactions between a 
fullercne molecule and one or two sodium atoms. The B3LYP hybrid method proposed by Becke 33 ~ 35 and included 
in the GAUSSIAN98 package 36 was used. This method is a mixture of Hartrec-Fock and DFT exchange terms 
associated with the gradient corrected correlation functional of Lee et al. 37 Calculations were made using LANL2DZ 
basis sets proposed by Hay et al. 38 where the double-zeta quality was used for the carbon all-electron calculations 
while effective core potentials were used for sodium atoms to replace the ten inner electron cores. For these atoms 
the valence electrons were described with double zeta quality basis sets proposed by the same authors. 

The first part of our calculation focussed on the adsorption of a single exohedral Na atom on the Ceo molecule. The 
DFT-LDA calculations performed by Hamamoto and coworkers 16 showed that the geometry of Cgo does not change 
significantly during the full geometrical optimization of Na„Cgo (with 1 < n < 12), compared to the case of pure 
Ceo- So, in the present study, the fullercne was frozen at its optimized geometry and, for NaCgo and Na2C@o, only 
the positions of the alkali atoms were relaxed. We can distinguish two kinds of CC bond lengths in the C 60 molecule. 
The simple C-C and the double C=C bonds used in this study were respectively frozen at 1.45 A and 1.37 A, close 
to the experimental values of Hedberg et al. 39 

B. Adsorption of a single alkali atom on Ceo 

We label H the adsorption site located above the center of an hexagonal site, P the site above a pentagon, BHH 
above the double C=C bond, BHP above the single C-C bond, and finally T above one carbon atom (vertex). In each 
case, only the radial distance d from the alkali atom to the centre of Ceo was optimized. The adsorption energies (E a ) 
of p sodium atoms on Cgo were calculated from the equation: 

£ a =£(C 60 )+px£(Na)-£(Na p C 60 ). (1) 

Different spin multiplicities were checked in order to get the most stable ground state. 

The adsorption energy, equilibrium distance, charge transferred and dipolc moments are presented in Table I for 
the five adsorption sites. As can be seen, H and P are the most stable configurations with a small preference of only 
0.03 eV for the hexagon. 

TABLE I: Adsorption energies E a , equilibrium distances Na-C, Mulliken charges on the alkali atom, and dipole moment of 
NaC6o- The present values are given in bold face. 



Adsorption site E a (cV) Na-C (A) q/e /j, (D) 

H 0.65 / 0.10°2.10* 2.74 / 5.08 a 2.69* 0.87 / 1.07* 14.53 / 14.63* 

P 0.62 / 0.10 a 2.05* 2.70 / 5.18 a 2.70* 0.87 / 1.06* 15.60 / 15.80* 

BHH 0.58 / 0.10 a 1.93* 2.59 / 5.00 a 2.47* 0.87 / 1.05* 17.08 / 16.54* 

BHP 0.55 / 0.10 a 1.78* 2.53 / 5.00 a 2.46* 0.85 / 1.06* 16.32 / 16.48* 

T 0.52 / 0. 10 a l. 98* 2.47 / 4.50 a 2.39* 0.85 / 1.04* 16.98 / 17.10* 

16.3±1.6 C 



"Unrestricted Hartree-Fock calculations [8] 
*DFT-LDA calculations [16] 
Experimental data [29,40] 

These results are in agreement with the calculations performed by Hamamoto et al. 16 who found H as the most 
stable configuration with a weak difference between H and P of about 0.05 eV. However, the adsorption energies 
reported in Ref. 16 are between 1.23 and 1.45 cV higher than ours. Since in both calculations double zeta basis 
sets were used, this difference is certainly due to the use of the local density approximation in the latter work. This 
approximation is known to often strongly overestimate adsorption energies. In contrast, Hira and Ray 8 performed 
unrestricted Hartree-Fock calculations of NaCgo and predicted a very weak adsorption energy of only 0.10 eV with 
no energy difference between the five different sites. Bond sites configurations, which are probably not true minima, 
show a small preference of 0.03 eV for the C=C bond, in agreement with the findings by Hamamoto and coworkers. 16 
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The less stable of the 5 configurations corresponds to the vertex site, indicating that the adsorption of a single alkali 
atom is favored by a maximal coordination. 

Estimates of the barrier for diffusion of Na over Cgo can be made from the data in Table I. We find barriers of 
0.07 eV between H sites or between P sites, and 0.07 and 0.10 eV between H and P sites. Thus we can expect a 
significant mobility at room temperature, and a possible influence of temperature on the average location of sodium 
atoms. 

The net charge on the alkali atom, as estimated from a Mulliken population analysis, 42 is also given in Table I. Our 
results are consistent with experimental 18 ' 28,40 and other theoretical 8,16 studies, which predicted an electronic transfer 
of about one electron from the singly occupied atomic orbital of the alkali atom to the lowest-unoccupied molecular 
orbital (LUMO) of C 60 . This behavior indicates a significant ionic bonding of the alkali atom with the fullerene. In 
the unrestricted Hartree-Fock calculations performed by Hira and Ray, 8 no significant electronic transfer between Na 
and Ceo was observed. This is correlated with a large distance (larger than 4.5 A) between Na and Ceo, and no actual 
bonding. 

The optimized bond lengths between the alkali atom and the nearest carbon atom calculated using the nonlocal 
functional are listed in Table I. These values show a decrease from H to P, then BHH and BHP, and finally to T sites. 
Thus the Na-C distance somewhat also reflects Na coordination. 

Finally, the dipole moments are given in Table 1 for all five configurations. Their values range from 14.53 D for site 
H to 17.08 D for site BHH. They agree well with the experimental data by Antoine et al., 40 who reported a dipole 
moment of 16.3 ± 1.6 D. Our results are also in good agreement with the DFT-LDA calculations of Hamamoto and 
coworkers. 16 



C. Adsorption of two alkali atoms on C 6 o 

We are now interested in the relative position of two alkali atoms adsorbed on the Cqo molecule. From results of 
the previous part (strong electronic transfer from alkali atom to Ceo molecule and weak diffusion barrier on the Cgo 
surface), we can expect that two alkali atoms will undergo a strong electrostatic repulsion. Two configurations were 
selected, with two adjacent sites or two opposite sites, respectively. Since the H sites are the most stable in the case 
of a single atom, we have restricted the present study to hexagonal sites only. In the case of opposite sites, only the 
radial distance to Cgo was optimized. For adjacent sites, the Na-Na distance was also optimized. 



TABLE II: Adsorption energies E a , equilibrium distances Na-C and Na-Na, Mulliken charges on each alkali atom, and dipole 
moment of Na2C6o- 



Adsorption site 


E a (eV) 


Na-C (A) 


Na-Na (A) 


q/e 


M(D) 


adjacent 


0.98 


2.55 


4.46 


0.83 


23.75 


opposite 


1.31 / 0.24 a 4.09 6 


2.72 / 4.48" 




0.85 / 1.05 6 






"Unrestricted Hartree-Fock calculations [9] 
''DFT-LDA calculations [16] 



The adsorption energies for Na2Ceo [Eq. (1)] are given in Table II. For the two configurations, adsorption energies are 
given along with the optimized Na-C and Na-Na distances, Mulliken charges on each alkali atom, and the associated 
dipole moment. The most stable configuration turns out to be the one with opposite sites, lower by about 0.34 eV 
than the other one. Even though the Mulliken charges on sodium are slightly smaller, they remain comparable to the 
single atom case. This significant electronic transfer (about 0.85e per atom) causes a strong electrostatic repulsive 
interaction between Na atoms, which in turn destabilizes the adjacent configuration. The present result is consistent 
with the DFT-LDA calculations of Hamamoto et al., 16 and with the experimental interpretations by Martin et al. 18 
for clusters with fewer than 7 sodium atoms. 



III. EMPIRICAL MODELLING 

The calculations presented in the previous section will now be used to parameterize an atomistic model, in order 
to achieve large scale sampling of the configuration space of larger clusters and perform simulations, which are even 
today beyond the possibilities of ab initio approaches. 
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A. Fluctuating charges potential 



The crucial element determining the relative extents of ionic-covalent forces and metallic forces is charge transfer 
between carbon and sodium. This transfer is expected to depend on how the metal atoms are located over Cgo- For 
instance, atoms closest to the C 60 surface are likely to undergo a stronger charge transfer than the most distant ones. 
Therefore we cannot use a model with fixed partial charges. A convenient answer to this problem is provided by the 
fluctuating charges model, often denoted as "fluc-q" potential. Originating from DFT analyses, 43 these ideas have 
been simultaneously but independently applied to a variety of diatomic molecules and some peptides by Rappe and 
Goddard, 44 as well as multicharged metal clusters by Sawada and Sugano. 45 More recently, the group of Berne 46 
has suggested how to improve molecular dynamics simulations with fluc-q potentials by using extended Lagrangian 
techniques, in a fashion similar to the Car-Parrinello method. 47 More rigorous derivations from density functional 
theory have since been proposed, 48,49 and also a treatment within hybrid quantum mechanical/molecular mechanical 
models. 50 The latest developments include higher-order electrostatic terms such as dipoles, 51 and the treatment of 
polarization forces. 52 Fluc-q potentials have been used by Ribeiro to study the dynamics in the glass-forming liquid 
Cao.4K .6(N03)i.4 53 as well as the Chcmla effect in molten alkali nitrates. 54 

Let us denote by R = {r^r^} the geometry of the Na„C 60 system, the vectors and representing the 3- 
dimensional positions of the n sodium atoms and 60 carbon atoms, respectively. We assume that the whole cluster 
carries a total charge Q, which may not necessarily be zero. The total potential energy V of the system is written as 
the sum of several contributions: 



V(R) = FNa({rJ) + ^ c ({r'}) + ^ nte r(R) + ^Coui(R). 



(2) 



Vn s in the pure metallic binding energy, which we took as an empirical many-body potential in the second moment 
approximation (SMA) to the electronic density of states in the tight-binding theory. 55 This potential involves five 
parameters £ 0j £ o, P, Q, and r a : 



2q(r ik /r -l) 



1/2 



(3) 



In this equation, the sum runs over all sodium atoms, and is the usual distance between atoms i and k. The initial 
values for the five parameters were taken from Ref. 56. Vq is the pure covalent binding energy of the Ceo molecule. 
This term was modelled with the Tersoff potential. 57 We did not attempt to change any of the original parameters, but 
we added several fixed point charges to improve the electrostatic description of this molecule. Following Schclkachclva 
and Tareyeva, 58 60 charges +Sq were placed on each C atom, while 30 charges — 25q were placed at the center of each 
C=C bond. Actually the full Tersoff potential is relatively costly, and we only used it for reoptimizing some stable 
structures. In the next section, the effects of approximating Cgo as rigid will be quantitatively addressed. 

The term T^ n t cr represents the non-ionic part of the interaction between carbon and sodium atoms. It is simply 
approximated as a pairwise repulsion, involving 2 new parameters D and (3: 



I'intor (R) 



E *>* 



iGNa„ jGCeo 



(4) 



We will assume in the following that chemical bonding between Na and C is controlled by charge transfer, hence we 
neglect other effects such as dispersion forces. In the fluc-q model, the Coulomb electrostatic interaction is expressed 
as 



F Coul (R) = E 



£Na<7i + 2^Na<7i 



+E 



E Jiilrfj 



E Jii'QiH' 



E Jjj'QjQj" 



t,3 



E 



(5) 



where and g' denote the charges carried by sodium atom i and carbon atom j, respectively. In Eq. (5) the labels 
and ) refer to sodium and carbon, respectively. £Na and eq are the electronegativities, i?Na and Hq the 
hardnesses, and the Coulomb interactions, which are assumed to be explicit in the interatomic distance ry. In 
practice, we have chosen the Ohno representation: 59 



J ij( r ) = V 2 + H ij 2 exp(-7ij' 



-1/2 



(6) 
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For a given geometry R, the effective electronegativities of each atom are defined by e = dV/dq, which yields 

£i = £Na + #Na<7i + ^ JijQ'j + ^ Jii'Qi' (?) 

for sodium atoms, and a similar equation for carbon atoms. Following Sanderson, 60 we now use the principle of 
electronegativity equalization. This principle states that, at equilibrium, all e's are the same in the molecule. If we 
denote by A this common value, then the charges {q i} q'j} are such that they minimize the potential Vc ou i above, under 
the constraint that the total charge is prescribed to a specific value Q. This minimization is equivalent to adding the 
extra term AQ^ 4 Qi + 1j ~ Q) t° the potential Vc ou i> and thus A plays the role of a Lagrange multiplier. 

Finding the optimal charges amounts to solving a (n + 61) x (n + 61) linear system. This system has a unique 
solution, as long as the hardnesses are all strictly positive. Otherwise the quadratic form Vcoui could no longer be 
minimized. 

We did not attempt to refine the present model further by adding higher order electrostatic or polarization terms 
and we assume that the self-consistency of charges is sufficient to ensure a correct electrostatic balance. Dispersion 
forces have been neglected as well. The interested reader is referred to the works in Rcfs. 51 and 52 for more 
details. However, we should mention that the fixed charges (8q, —2Sq), denoted as qk in what follows, contribute to 
a better electrostatic description for the whole system. The fluc-q potential can easily account for such extra fixed 
charges as well as for a possible external electric field E. Both fixed charges and electric field affect the self-consistent 
determination of the fluctuating charges, and the potential Vboui should be supplemented with the term 

JikQilk + Jj,kq'j<lk + Y ( E ' r ^ qi + X!( E ' r ^ q r ( 8 ) 

i,k j,k i j 



B. Parameterization 



The empirical fluc-q potential described in the previous paragraph has 15 parameters, namely the five parameters 
of the SMA potential, the two parameters of the Na-C repulsion, the electronegativity difference ec — £Na, the six 
parameters jij and Hij, which include the two hardnesses i?Na and Hq, and finally the fixed charge 5q. These 
parameters were chosen so as to reproduce several properties previously computed by DFT calculations, or measured 
in experiments. However these data do not give much insight into the bonding between Na and C atoms. In addition, 
some of the parameters, especially those involved in the Jij(r) functions, can be in principle extracted from calculations 
on diatomic molecules only. Also, while Na2 and C2 are well characterized experimentally and theoretically, almost 
no data is available for the NaC molecule. 

We have therefore achieved ab initio Configuration Interaction (CI) calculations of NaC, as well as some calculations 
on Na2 and C2 to get the Coulomb integrals. These were taken in minimal linear combinations of atomic orbitals basis 
set obtained by contracting the Gaussian-type orbitals of atomic self-consistent field calculations in the previous basis 
sets. The pointwise integrals were then fitted through the Ohno functions, Eq. (6) above. The ab initio calculations 
on Na2 also provided estimates of the binding energy and equilibrium distance that should be predicted by the model. 



TABLE III: Data used for fitting the empirical model with fluctuating charges. 



Property 



Source 



Reference 



Predicted 



q/e (NaCeo) 
fj, (NaCeo, D) 
AE (Na 2 C 6 o, eV) 
r(Na 2 , A) 
E(N& 2 , eV) 
a(C 60 , A 3 ) 



B3LYP 
B3LYP 
B3LYP 

ab initio" 
ab initio" 
exp. 6 



0.87 
14.5 
-0.34 
3.06 
0.74 
76.5 



"Data from Jeung, Ref. 61 . 

'Data by Antoine et al. from Ref. 62 



0.88 
14.5 
-0.27 
3.34 
0.58 
58.5 



Appropriate fits of the Coulomb integrals gave us initial guesses for the parameters 7^ and . We then obtained 
the whole set of parameters by minimizing a standard error function x 2 , to reproduce the following quantities, quoted 
in Table III: 



6 



1. the charge transfer from sodium (0.87e) and the electric dipole (14.5 D) in NaCgo, values taken from our DFT 
calculations; 

2. the energy difference between the two Na2C6o isomers with sodium atoms on adjacent or opposite hexagonal 
sites (AE = —0.34 eV from our DFT calculations, in favor of the opposite sites location); 

3. the binding energy (0.74 eV) and equilibrium distance (3.06 A) in Na2 from ab initio calculations; 61 

4. the experimental electric polarizability of C 60 ( 76.5 A 3 , taken from Ref. 62). 

The error function \ 2 was then minimized using a local Monte Carlo search. The moves in the parameters space 
were not taken too large, in order to keep the parameters close to their initial magnitudes. The final values of the 
parameters of the SMA potential are e = 0.01366 eV, £ = 0.1689 eV, p = 7.175, q = 1.34, and r = 3.29 A. The 
parameters for the repulsive potential between Na and C atoms are D = 313.4 eV and (3 = 3.423 A -1 . For the 
Coulomb integrals, we found # Na = 7.16 eV, H c = 13.68 eV, # Na c = 6.32 eV, 7 Na = 0.129 A~ 2 , 7c = 1.02 A~ 2 , 
and 7NaC = 0.15 A~ 2 . Finally, the electronegativity difference was taken as Ae = 4.6 eV, and the fixed charge Sq as 
0.209. We represented in Fig. 1 the Coulomb interactions between Na and C atoms. They are found to remain very 
close to the ab initio calulations. 




0.0 1 ' 1 ■ 1 ' 1 ' 1 ' 1 

2 4 6 8 10 

r (bohr) 

FIG. 1: Coulomb integrals for the Na2, C2 and NaC molecules. The curves are given for the empirical potential (thick lines) 
and for the ab initio calculation (thin lines). 

The predicted values of the quantities used in the fitting procedure are also reported in Table III. The general 
consistency is very good, even though we note that the polarizability of Ceo is slightly underestimated in our model. 
However, only the Sq parameter has some influence on this parameter. The energy difference between the two isomers 
of Na2C6o is somewhat overestimated, but the crucial point here is that the isomers are ordered similarly. 

C. Simulation tools 

The putative global minima, or lowest energy structures, were determined using a variant of the Monte 
Carlo+minimization algorithm, 63 also known as basin-hopping. 63,64 This optimization stage was made with a fully 
rigid Ceo molecule, and the random atomic displacements of the sodium atoms were of two kinds. With a probability 
p, the selected atom is moved from its equilibrium position by a vector Ar with maximum amplitude 10 A. With 
probability 1 — p, it is moved anywhere over the Ceo with the same radial distance from the centre of mass of the 
fullerenc. In practice, all sodium atoms are moved before a quench is performed to find a new stable minimum. The 
value of p was taken as 80% for the smaller sizes, and decreased down to 10% for n > 12. To prevent divergences 
in the fluc-q potential due to short distances resulting from the random displacements, and to prevent pathological 
crossing of the Ceo surface by Na atoms, we added a purely repulsive potential between the surface of Ceo and the 
sodium atoms: 

^cp({r J }) = ^|(||r J |t-r C60 ) 2 , (9) 
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where we took n — 880 eV/A 2 and rc 60 = 3.89 A. This potential was only used in concern with global optimization. 

We also carried out some finite temperature molecular dynamics simulations using the Nose- Hoover chain thermostat 
technique. 65 In the extended Lagrangian scheme, the charges were also thermostatted at a low temperature T* = 
T/100, T being the vibrational temperature, to prevent divergence of their kinetic temperature. 66 The masses of 
the thermostats were taken as Qi = 3nksT/Ld 2 (1st thermostat) and Qk>i = Qi/3n (other thermostats) for the 
atoms, and Q\ = (n + 60)fc_BT*/w 2 (1st thermostat) and Q*k>i = Q\/{n + 60) (other thermostats) for the charges. 65 
Appropriate values for the square frequencies were taken as J 2 = 10~ 6 and Wq=5x 10~ 6 atomic units, respectively. 
The fictitious masses of the charges were chosen equal for all charges, and taken as 10 3 au. The velocity Verlet 
algorithm was used to propagate the dynamics with a time step At = 1 fs, and the parallel tempering strategy 67 was 
used to accelerate convergence when simulating the clusters at thermal equilibrium. 

IV. RESULTS 
A. Structural properties 

For each size n in the range 1 < n < 30, 10 4 quenches were performed during the basin-hopping searches for the 
global minimum. The most stable structures found using this algorithm are reported in Fig. 2. We have also indicated 
the (non-Ci) symmetry groups. Although this may not be obvious from looking at Fig. 2, the Cgo molecule is also 
optimal in these structures when modelled with the Tersoff potential. 57 The effect of adding sodium atoms on the Ceo 
geometry will be quantified later. 

At low coverage, Na atoms tend to stay as far as possible from each other, in order to minimize Coulomb repulsion 
induced by a significant charge transfer. The growth roughly goes symmetric until the seventh atom is added. There 
are then only few open spaces for the eighth atom to stay, and the structure found by breaking some Na-C bonds in 
favor of metal Na-Na bonds appears slightly more stable. More importantly, this droplet initially created for n = 8 
remains at larger sizes and grows monotonically as further atoms are deposited. Actually the droplet even captures 
some of the isolated atoms during the growing process: only two atoms remain isolated for n = 30. 

These are the main results of the present work. Our model predicts that charge transfer is initially strong enough to 
inhibit metallic bonding. Then, as more atoms are added, the competition between ionic and metal forces finally turns 
at the advantage of the latter. The crossover size for this wetting-to-nonwetting' transition, which we estimate here 
to be around n = 8, may be of course somewhat different in more accurate chemical descriptions. It is however worth 
pointing out that, even though the present potential is not explicitely quantal, it still shows an enhanced stability 
of Na trimers upon nucleation. The atom unconnected to the Ceo carries a slightly negative charge, and the whole 
trimer nearly has the charge +1, in agreement with the known special stability of Nag". 68 

In Fig. 3 we have plotted two indicators of the mistake made when considering the Cgo molecule as rigid. The 
whole Na„C6o systems were locally reoptimized using the Tersoff potential, and we calculated the relative error 
e = L4(rigid) — A(Tersoff)]/A(Tersoff) for two properties A: (1) the global shape is estimated from A = ^ ||i"i|| 2 , 

i 

where the sum runs over all Na and C atoms, and (2) the total charge carried by sodium atoms, A = qi. 

As can be seen in Fig. 3, the relative error remains smaller than 5 x 10~ 4 for the shape, and smaller than 2 x 10~ 4 
for the charge. Of course, heating the cluster would introduce some degree of fioppyness in the fullerene. However, 
at the temperatures involved in the experiment (the order of magnitude of 300 K), we do not really have to worry 
about the vibrational motion of carbon atoms. The Cgo molecule itself only melts above 2000 K. 69 Hence it is quite 
safe to use the rigid approximation in the present work. We note, however, that other materials such as fluorine may 
spontaneously induce significant deformations when they cover Ceo in large amounts. 70 

B. Electrostatics 

The fluctuating charges potential can be straightforwardly used to compute electrostatic properties, some of them 
being amenable to experimental comparison. The total charge carried by all sodium atoms is represented in Fig. 4 
versus cluster size. As long as new hexagonal sites of Cgo are capped by the added Na atoms, charge transfer increases 
steadily. The total charge reaches a plateau at n — 8, corresponding to the onset of nucleation. Above this size, slightly 
more than 3 electrons are transferred in the present empirical model. Single point DFT calculations performed on 
some of the smaller sizes confirm this trend, but they estimate the maximal charge transfer to be closer to 5-6 rather 
than 3. Such a larger value is in agreement with the theoretical findings of Hamamoto and coworkers. 16 This may 
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(7,C 3 ) (8) (9) (10) (11) (12,C S ) 




(13) (14) (15) (16) (17) (18) 




(19) (20) (21) (22) (23,C S ) (24) 




(25) (26) (27) (28) (29) (30) 




FIG. 2: Lowest-energy structures found for Na„C60 in the range 1 < n < 30 using the fluctuating charges potential. 

be associated with the role in the molecular orbital calculations of the threefold degenerate LUMO orbital of Ceo, an 
effect which cannot be explicitely accounted for in the present continuous electrostatic model. 

In Fig. 4 we also plotted the variations of the magnitude of the electric dipole moment /i. Two different regimes 
are again observed, that can be correlated with the known structural transition. Below 8 sodium atoms, the rather 
symmetric structures due to Coulomb repulsion exhibit some odd-even alternance in the dipole. From n — 8 and 
beyond, the growing droplet shows a much more progressive increase of the dipole, except close to sizes that lose one 
of the isolated atoms. 

The zero temperature electric polarizability a was calculated by averaging the diagonal part of the corresponding 
3x3 tensor. Small electric fields (10 -6 au) were added along the three Cartesian axes, and the dipoles were obtained 
after reoptimization. At finite temperature corresponding to the experiments by Dugourd et al., 29 we estimated the 
electric susceptibility \ by assuming that the dipoles were rigid and statistically oriented within the electric field. 
This allows us to use a simple Langevin formula for the susceptibility, namely \ = a + M 2 /3fcsT. The variations of a 
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FIG. 3: Relative error for the geometric distortion (empty squares) and for the total charge on sodium (full circles) when 
considering the Ceo molecule as rigid, with reference to the Tersoff potential. 
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FIG. 4: Total charge (empty circles) and electric dipole moment (full squares) versus size in Na„C6o clusters. 



and x(300 K) with size are represented in Fig. 5 for the entire range 1 < n < 30. 

The polarizability shows no particular dependence with n besides a steady but slow increase. On the contrary, 
because /z has large variations with size, this also holds for \- Experimental data 29 are in semi-quantitative agreement 
with the present results. While the general behavior observed by Dugourd et al. 29 is similar to the one for \ in Fig. 5, 
the values they measure are significantly larger. This might be due to the neglect of probable increases in /j, as 
temperature increases, due to the floppyness of these clusters at 300 K. 



C. Charge and temperature effects 



In experiments, clusters are ionized for subsequent size-selection. To see the general influence of ionization on the 
growing process of Na clusters over Ceo, we have performed additional global optimization on anionic (Q = —1) and 
cationic (Q = +1) systems. The size of the largest sodium fragment is represented in Fig. 6 versus the total number of 
sodium atoms. Here a fragment is a set of connected atoms, where any two atoms are said to be connected whenever 
they are distant by less than 8 bohr. This quantity is suitable for detecting the onset of nucleation. Indeed, we can 
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FIG. 5: Electric polarizability (empty circles) and susceptibility at 300 K (full squares) versus size in Na n C60- 



easily identify in Fig. 6 the sizes of neutral clusters where the droplet captures one of the remaining isolated atoms. 
These sizes are found for n = 13, 18, and 22. 




Cationic clusters are less stable than neutrals, as they start to nucleate a droplet only above n — 9. This should be 
attributed to the larger coulombic repulsion the new atoms undergo as they are placed close to already present atoms. 
On the contrary, charging negatively the system diminishes charge transfer on sodium and therefore helps in creating 
metallic bonds. This is precisely what we see in Fig. 6, where the crossover size between coating and segregation is 
only 3. Obviously, the above arguments should be further discussed in the light of quantum effects. 

The finite temperature equilibrium properties have been simulated using constant temperature, extended Lagrangian 
molecular dynamics as described in the previous section. The simulations were carried out for 30 simultaneous 
temperatures equidistant by AT = 10 K, and were each 10 ns long, after 2 ns were initially discarded for equilibration. 
To prevent the alkali atoms from dissociating at high temperatures, we enclosed the whole system inside a spherical 
container centered around the Ceo- 

The time averaged size of the largest fragment is also represented in Fig. 6, at T = 300 K, for 1 < n < 30. This 
quantity essentially follows the K behavior, but is slightly higher than 1 for n = 6 and 7. For n > 8, the 300 K value 
is slightly lower than the static value. This suggests that some metallic bonds are weak and easily broken for larger 
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droplets. One must then recall that 300 K is quite a large temperature for sodium clusters, and that dissociation 
events would tend to occur on the droplet is there were no artificial container to prevent it. 

The influence of the fullerene of the global finite-temperature behavior of the sodium cluster was investigated 
by calculating the heat capacity of Na 2 oC6o, and comparing to our previous results on bare Na clusters. 71 The 
thermodynamical data was analysed using a multiple- histogram reweighting procedure. 72 The heat capacity curves 
are given in Fig. 7, and we also represented the variations of the (arithmetic) average fragment size (n) with increasing 
temperature. 

60 



50 

cq 

^ 40 
30 



20 

100 200 300 

T(K) 

FIG. 7: Heat capacities of Na2o in gas phase or deposited on Ceo- Inset: average fragment size (n) versus temperature in 
Na2oC6o- 

The heat capacity of Na2oC6o looks slightly different from that of Na2o, with an extra shoulder bear 40 K. The main 
features at ~ 100 K and ~ 200 K in the bare cluster are enhanced in the supported cluster. Looking at snapshots 
of configurations taken from the T — 50 K trajectory indicates that some sodium atoms can change adsorption sites, 
especially isolated ones. The peak at 100 K can be correlated with prcliminar isomcrization of the Nai7 droplet, from 
which some atoms can dissociate while skating over the C 6 o surface. This is best evidenced on the small decrease 
of (n) with T. On the contrary, at 180 K an inverse process takes place, where the previously isolated atoms grow 
onto the droplet. Additionally some Na 2 dimers may be spontaneously formed, and briefly destroyed, during the 
simulation. This further increases the average fragment size. 

D. Electric field 

The experiments performed by Dugourd and coworkers 29 involve the clusters travelling through a small but very 
intense region of active electric field. Actually, both the field and its gradient are strong. Because the C and Na atoms 
get charged differently when put into contact in Na„C6o they may react differently in presence of such a strong field. 
We attempted to quantify these effects by carrying out extra global optimization for Na 2 C6o within a field of constant 
magnitude E = 5 x 10~ 4 au, corresponding to 2.55 x 10 8 V m _1 . In addition to the 4 atomic degrees of freedom 
(neglecting the radial distance to the carbons), the direction of the field must also be optimized. Using this approach 
we found 3 isomers to be potentially the most stable when E is waried between and 10~ 3 au. These isomers are 
represented in Fig. 8 along with the variations of their energies with increasing field. In isomers B and especially C 
the two alkali atoms get closer to each other. This is not surprising, since the presence of a strong electric field causes 
the positively charged ions to move in its direction. As a matter of fact, the field is optimally perpendicular to Na 2 
in both B and C, while it has little effect when fully aligned (isomer A). Above E ~ 10~ 4 au, isomer B becomes the 
most stable until E reached 2.5 x 10~ 4 , above which C is the new global minimum. 

The present results therefore show that a strong, external electric field similar to the one used in experiments can 
indeed induce some configurational changes in the structure of Na„ adsorbed on Ceo- Of special interest to us, the 
field tends to favor metallic bonding by lowering the effects of ionic repulsion. Hence we can expect the onset of 
nucleation to be earlier in the presence of a field. 
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FIG. 8: Variations of the binding energy of three isomers Na2C6o versus the magnitude of electric field. The isomers (A, B, 
and C) are represented on the lower panel, the arrow indicates the optimal field direction. 

E. Dynamics 

For each temperature in the range 50 K< T <300 K with steps AT = 50 K, we have performed 5000 short (10 ps 
long) MD trajectories from which we computed the average atomic mean square displacement, and subsequently the 
diffusion constant D. This was repeated for a number of sizes, and also for a non-zero external electric field. The 
data for n = 1, 4, 6, 12, and 20 are represented as an Arrhenius plot in Fig. 9. The general linear variations of log!? 
versus 1/T are characteristic of diffusive motion. The activation barriers extracted from these plots are approximately 
A ~ 400 K at zero field, and A ~ 600 K at E — 5 x 10~ 4 au. At room temperature sodium atoms thus show a 
significant mobility, whatever a droplet is present or not. Adding an electric field acts as a quench of the dynamics 
which could trap the system into metastable configurations, provided that the vibrational modes of Ceo take a part 
of the excess energy. 

V. DISCUSSION AND CONCLUSION 

In this work, we have performed ab initio and DFT calculations on NaC, NaCgo and Na2Cgo in order to construct 
en empirical model for larger clusters. This model allowed us to achieve unconstrained global optimization as well as 
(thcrmo)dynamical investigations in a broad size range. Our main findings are as follows: 

1. Based on first-principles calculations, Na2 is most stable in dissociated form, the two sodium atoms lying on 
opposite hexagonal sites of the fullerene; 

2. In its initial stages, the growth of Na„ on Cgo occurs by minimizing Coulomb repulsion between the nearly 
cationic alkali atoms, thus placing them as far as possible from each other; 

3. At n ~ 8, the electrostatic penalty is overcome by creating metallic bonds, and a sodium droplet is formed. The 
droplet grows for larger sizes; 
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FIG. 9: Diffusion constants of various Na„C60 clusters versus inverse temperature. For NaeCeo, the results of simulations with 
non-zero external electric field are also plotted for comparison. 



4. Depending on the ionic nature of the cluster, the crossover between homogeneous coating and droplet formation 
varies from 3 (anion) to 10 (cation); 

5. A finite temperature or a finite external electric field can enhance metallic bonding and further decrease the 
crossover size; 

6. At room temperature, the alkali atoms undergo a significant mobility over the C@o surface. 

From these results, a two-stage coating process emerges. This picture is consistent with the observations by the 
Martin group 18 and by Palpant et al. 26 ' 28 who interpreted their mass spectrometry and photoclcctron spectroscopy 
data as the signature of homogeneous coating at low sizes. However, in accordance with the electric susceptibility 
measurements by Dugourd and coworkers, 2 we find that larger clusters are more stable when they form a droplet. 
Actually, we predict that Na n C6o clusters with 10 < n < 30 show an intermediate structure, which consists of a main 
droplet and some remaining isolated atoms. 

Even though the empirical model was parameterized on ah initio calculations and experimental data on small 
clusters, it is nothing more than a model, and should thus be considered as an approximation, especially when used 
for larger clusters. Therefore, it may well be that the transition between wetted and non-wetted forms takes place at 
slightly different size if we include other effects such as polarization, or if we describe metallic bonding using a more 
realistic Hamiltonian such as tight-binding. Also, the existence of rather large sodium coverage of opposite isolated 
atoms, which mainly results from electrostatic balance in avoiding too large dipole moments, might disappear if more 
realistic screening due to polarization forces or quantum effects is involved. This is to be further checked. 

Nevertheless our conclusions help in solving the experimental puzzling results of various authors. Furthermore, 
it is remarkable that very recent photodissociation and photoionization experiments carried out by Pellarin and co- 
workers 81 on the same systems reached the same qualitative and quantitative conclusions as ours. At low coverage 
(n < 6) they observe ionic-like bonding and strong charge transfer. As the number of metal atoms increases, metallic 
bonds appear and become preponderant, finally controlling cluster properties. The picture of a metallic cluster 
deposited on the Cqo molecule is supported by the odd-even alternation seen in the stability pattern, where even 
(resp. odd) clusters evaporate mostly single atoms (resp. dimers). 

The latter effect is essentially quantum mechanical and lies beyond the empirical potential described in the present 
work. The favored formation of trimers at low coverage, suggested by Palpant and coworkers 28 and further supported 
by the calculations of Hamamoto et al., 16 is not reproduced here. This is partially due to the neglect of quantum 
effects, but also to the strong electrostatic repulsion between the alkali atoms. At the onset of nucleation, a trimer 
is actually formed, which lies perpendicular to the C@o surface. Examination of the partial charges shows that, while 
the two atoms in contact with Ceo are significantly positively charged, the most outer atom is charged negatively, 
therefore enhancing the stability of the trimer. It would be interesting to check whether this effect is also seen in 
electronic structure calculations. 

Finally, we would like to stress that the present model is not limited to simple metals adsorbed on fullerenes. 
Upon careful parameterization using experimental or first-principles data, it could be suitably modified to treat any 
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other metal, provided that the specific metallic interaction is taken into account. It could also be used for other 
materials, which are of experimental or technological interest. For instance, the single-crystal X-ray structure of 
C60F48 determined by Troyanov and coworkers 70 could be reproduced using a similar bond-order potential including 
fluctuating charges by Stuart et al. 73 More generally, it is appropriate for treating complex systems, which involve 
partially covalent and metallic bonding where ionic effects and charge transfer can prove to be a determining factor. 
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APPENDIX: DETAILS OF THE CONFIGURATION INTERACTION CALCULATION FOR NaC 

Some information about the distance dependence of the interaction between Na and C in NaCeo can be extracted 
from the NaC molecule. This must however be done with care, since charge transfer is intricated with the multiplct 
structure of carbon. In order to get a global picture, we have thus determined not only the ground state but also the 
lowest valence excited states. 

The calculation was done representing both C and Na via standard valence semi-local pseudopotentials of the 
Barthelat and Durand type 74 with [Is 2 ] and [ls 2 2s 2 2p 6 ] cores respectively. The calculation was carried out within 
the Linar Combination of Atomic Orbitals (LCAO) expansion scheme with uncontracted gaussian type functions, 
namely 5s/5p/5d on both centers in order to decribe the valence electrons. A core polarization operator of the form 
— l/2a c / c / c was added on the sodium core in order to take into account the core- valence polarization and correlation 
following to the formulation of Miiller and Meyer 77 for the electronic contribution to the electric field / on sodium. 
The Na + core polarizability was a c = 0.995 a 3 ,, and a unique stepwise cut-off radius 76 of 1.45 ao was used. The 
core-polarization contribution on carbon was neglected due to the extremely small value of the dipole polarizability 
of the C 4+ ion with respect to that of Na + . 

The CI treatment was achieved using the multi-reference variational-perturbative CI algorithm CIPSI 80 within its 
Quasi-Degenerate Perturbation Theory version 79 and the Mollcr-Plcsset partition of the Hamiltonian. The intermulti- 
plet separations of the carbon atom, calculated consistently, were found to be A( 3 P— 1 D) = 1.386 eV, A( 3 P— 1 5)=2.565 
eV, versus 1.260 and 2.680 eV experimentally. 78 For C~ the separations were found to be A( 4 P — 2 £>)=1.656 eV, 
A( 4 P- 2 5)^2.095 eV. 

The lowest states potential energy curves are shown in Fig. 10. Combination of Na(3s 2 S) with C( 3 P) generates 
three molecular states, namely 1 4 £~, 1 2 £~ and 1 2 II. Combination with C( 1 £>) also generates three states, namely 
1 2 A~, 2 2 II and 1 2 S + . Finally, combination with C( 1 S I ) generates a single state 2 2 S + . In the following, the valence 
molecular orbitals (MO's) will be labelled lc, 2<r, 3a, lir x , and ln y . Despite hybridization, they can be qualitatively 
associated with the 2s, 2p z , 3s, 2p x and 2p y atomic orbitals, respectively. 

The ground molecular state correlated with the C( 3 P) asymptote is a quartet state, 1 4 E~. At its equilibrium 
distance i? e =4.40 ao, it is mainly spanned by molecular configuration \a 2 2u 1 \'K 1 x 'K v with all the electrons in the 2p 
shell of carbon and is therefore expected to be strongly ionic. The DFT-B3LYP calculation of the NaC ground state 
provides very similar results, see Fig. 10. In particular, the dipole moment is found to be 7.88 D at the equilibrium 
distance. This is confirmed by the dipole moment value /x = 8.53 D at the equilibrium distance and the associated 
charge transfer 5q = 0.80. It can be clearly seen in Fig. 10 that the ground state adiabatically undergoes at R = 7.6 ao 
a strong avoided crossing with a state of same symmetry 2 4 £~. As a consequence of the configuration switch, the 
adiabatic state 2 4 £~ has a small dipole moment in the short distance range R < 6 ao. In the long distance range 
10 ao < R < 14 ao, it is correlated with C~( 4 P) + Na + and its potential curve exhibits a —1/R behavior, consistent 
with a continuously increasing dipole moment (up to \i « —21.6 D at R = 12 ao, the value for a perfectly ionic state 
would be —30.5 D). For still larger distances, the state interacts with another covalent excited state. The strong 
coupling between the two quartet states 1 4 £~ and 2 4 £~ states explains the significant dissociation energy of the 
ground state D e = 1.464 eV. 

Most other states dissociating into neutral asymptotic limits show similar (sometimes multiple) avoided crossing, 
mainly due to their interaction with ionic configurations. The positions of the avoided crossings are obviously deter- 
mined by the respective positions of the multiplcts of carbon and those of the ionic states dissociating into C~+Na. 
This also influences the dissociation energies of the lowest states stabilized by the interactions. In consistency with this 
picture, one can notice the avoided crossing between states 1 2 £~ and 2 2 S— , those between the three 2 II states and 
finally that between 1 2 A and 2 2 A. As a result, the lowest states involved in the avoided crossings are stabilized. This 
leads to significant dipole moments and charge transfer values at equilibrium for two other states correlated with the 
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FIG. 10: Potential energy curves of the NaC molecule. All continuous curves correspond to CI calculations. Some DFT values 
are also reported (dots) for the ground state. 



lowest asymptote, and explains their bonding properties (D e = 0.30 eV, R e — 4.6 a for state 1 2 £~, D e = 0.60 eV, 
R e = 4.7 a for state 1 2 IT). The same holds for state 1 2 A (D e = 1.61 eV, R e = 4.4 a ) dissociating into C( 1 D). State 
1 4 II is the only exception to the previous behavior, since it is spanned by configuration lo' 2 2<7 1 l7r;J,l7r 2 3(7 1 . Indeed 
within a single configuration approximation, no charge transfer can take place between the 2a and 3ct MO's in the 
quartet state because all singly occupied orbital spins are aligned. This state can thus be described as essentially 
covalent with a weak dipolc moment and small charge tranfer (Sq = 0.25 at R = 4 ao). Its potential curve is essentially 
non bonding and can be fitted via a simple exponential form, which can be used as the NaC covalent repulsive term 
of the model, Eq. (4). 
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